Fit VBGE models to back-calculated length-at-age

Load libraries

pkgs <- c("tidyverse", "tidylog", "rTPC", "nls.multstart", "broom", "RColorBrewer", 
          "viridis", "minpack.lm", "patchwork", "ggtext", "brms", "tidybayes", "modelr")

# minpack.lm needed if using nlsLM()
if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library,character.only = T))

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'tidylog'


The following objects are masked from 'package:dplyr':

    add_count, add_tally, anti_join, count, distinct, distinct_all,
    distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
    full_join, group_by, group_by_all, group_by_at, group_by_if,
    inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
    relocate, rename, rename_all, rename_at, rename_if, rename_with,
    right_join, sample_frac, sample_n, select, select_all, select_at,
    select_if, semi_join, slice, slice_head, slice_max, slice_min,
    slice_sample, slice_tail, summarise, summarise_all, summarise_at,
    summarise_if, summarize, summarize_all, summarize_at, summarize_if,
    tally, top_frac, top_n, transmute, transmute_all, transmute_at,
    transmute_if, ungroup


The following objects are masked from 'package:tidyr':

    drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
    spread, uncount


The following object is masked from 'package:stats':

    filter


Loading required package: viridisLite

Loading required package: Rcpp

Loading 'brms' package (version 2.19.0). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').


Attaching package: 'brms'


The following object is masked from 'package:stats':

    ar



Attaching package: 'tidybayes'


The following objects are masked from 'package:brms':

    dstudent_t, pstudent_t, qstudent_t, rstudent_t



Attaching package: 'modelr'


The following object is masked from 'package:broom':

    bootstrap

# devtools::install_github("seananderson/ggsidekick") ## not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Load functions
home <- here::here()
fxn <- list.files(paste0(home, "/R/functions"))
invisible(sapply(FUN = source, paste0(home, "/R/functions/", fxn)))

Load cache

# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/analyze-data/02-fit-vbge_cache/html"))

Read and trim data

d <- #read_csv(paste0(home, "/data/for-analysis/dat.csv")) |> 
  read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv") |> 
  filter(age_ring == "Y", # use only length-at-age by filtering on age_ring
         !area == "VN") 

Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
filter: removed 36,279 rows (11%), 302,181 rows remaining

# Sample size by area and cohort
ns <- d |> 
  group_by(cohort, area) |> 
  summarise(n = n())

group_by: 2 grouping variables (cohort, area)
summarise: now 453 rows and 3 columns, one group variable remaining (cohort)

# Minimum number of observations per area and cohort
d$area_cohort <- as.factor(paste(d$area, d$cohort))

d <- d |>
  group_by(area_cohort) |> 
  filter(n() > 100)

group_by: one grouping variable (area_cohort)
filter (grouped): removed 2,796 rows (1%), 299,385 rows remaining

# Minimum number of observations per area, cohort, and age
d$area_cohort_age <- as.factor(paste(d$area, d$cohort, d$age_bc))

d <- d |>
  group_by(area_cohort_age) |> 
  filter(n() > 10)

group_by: one grouping variable (area_cohort_age)
filter (grouped): removed 3,584 rows (1%), 295,801 rows remaining

# Minimum number of cohorts in a given area
cnt <- d |>
  group_by(area) |>
  summarise(n = n_distinct(cohort)) |>
  filter(n >= 10)

group_by: one grouping variable (area)
summarise: now 11 rows and 2 columns, ungrouped
filter: no rows removed

d <- d[d$area %in% cnt$area, ]

# Plot cleaned data
ggplot(d, aes(age_bc, length_mm)) +
  geom_jitter(size = 0.1, height = 0, alpha = 0.1) +
  scale_x_continuous(breaks = seq(20)) +
  theme(axis.text.x = element_text(angle = 0)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  labs(x = "Age", y = "Length (mm)") +
  guides(color = "none") + 
  facet_wrap(~area, scale = "free_x")

# Longitude and latitude attributes for each area
area <- c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA", "TH")
nareas <- length(area)
lat <- c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4, 56.1)
lon <- c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7, 15.9)
area_attr <- data.frame(cbind(area = area, lat = lat, lon = lon)) |>
  mutate_at(c("lat", "lon"), as.numeric)

mutate_at: converted 'lat' from character to double (0 new NA)
           converted 'lon' from character to double (0 new NA)

Fit VBGE models

# Get individual growth parameters (functions: VBGF/Gompertz,nls_out,fit_nls)
IVBG <- d |> 
  group_by(ID) |> 
  summarize(k = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$k,
            k_se = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$k_se,
            linf = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$linf,
            linf_se = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$linf_se)

group_by: one grouping variable (ID)

summarize: now 79,122 rows and 5 columns, ungrouped

Inspect predictions

IVBG <- IVBG |> drop_na(k) # The NAs are because the model didn't converge or because they were below the threshold age

drop_na: removed 41,442 rows (52%), 37,680 rows remaining

IVBG <- IVBG |>
  mutate(k_lwr = k - 1.96*k_se,
         k_upr = k + 1.96*k_se,
         linf_lwr = linf - 1.96*linf_se,
         linf_upr = linf + 1.96*linf_se,
         row_id = row_number())

mutate: new variable 'k_lwr' (double) with 37,641 unique values and 0% NA

        new variable 'k_upr' (double) with 37,641 unique values and 0% NA

        new variable 'linf_lwr' (double) with 37,641 unique values and 0% NA

        new variable 'linf_upr' (double) with 37,641 unique values and 0% NA

        new variable 'row_id' (integer) with 37,680 unique values and 0% NA

# Plot all K's
IVBG |>
  #filter(row_id < 5000) |>
  ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr)) +
  geom_point(alpha = 0.5) +
  geom_errorbar(alpha = 0.5) +
  NULL

# Plot all L_inf's
IVBG |>
  #filter(row_id < 5000) |>
  ggplot(aes(row_id, linf, ymin = linf_lwr, ymax = linf_upr)) +
  geom_point(alpha = 0.5) +
  geom_errorbar(alpha = 0.5) +
  NULL

# We can also consider removing individuals where the SE of k is larger than the fit
IVBG |>
  #mutate(keep = ifelse(k > quantile(k_se, probs = 0.95), "N", "Y")) |>
  mutate(keep = ifelse(k < k_se, "N", "Y")) |>
  #filter(row_id < 10000) |>
  ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr, color = keep)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~keep, ncol = 1) +
  geom_errorbar(alpha = 0.5) +
  NULL

mutate: new variable 'keep' (character) with 2 unique values and 0% NA

# Add back cohort and area variables
IVBG <- IVBG |> 
  left_join(d |> select(ID, area_cohort)) |> 
  separate(area_cohort, into = c("area", "cohort"), remove = TRUE, sep = " ") |> 
  mutate(cohort = as.numeric(cohort))

Adding missing grouping variables: `area_cohort_age`

select: dropped 11 variables (length_mm, age_bc, age_catch, age_ring, area, …)

Joining with `by = join_by(ID)`
left_join: added 2 columns (area_cohort_age, area_cohort)
> rows only in x 0
> rows only in y ( 97,478)
> matched rows 198,323 (includes duplicates)
> =========
> rows total 198,323
mutate: converted 'cohort' from character to double (0 new NA)

# Summarise and save for sample size
sample_size <- IVBG |>
  group_by(area) |> 
  summarise(n_cohort = length(unique(cohort)),
            min_cohort = min(cohort),
            max_cohort = min(cohort),
            n_individuals = length(unique(ID)),
            n_data_points = n())

group_by: one grouping variable (area)
summarise: now 11 rows and 6 columns, ungrouped

sample_size

# A tibble: 11 × 6
   area  n_cohort min_cohort max_cohort n_individuals n_data_points
   <chr>    <int>      <dbl>      <dbl>         <int>         <int>
 1 BS          18       1985       1985          2370         11832
 2 BT          27       1972       1972          1882          9055
 3 FB          48       1969       1969          8678         47933
 4 FM          38       1963       1963          6704         39230
 5 HO          34       1982       1982          2711         12454
 6 JM          63       1953       1953          7738         40229
 7 MU          19       1981       1981          1920          9930
 8 RA          15       1985       1985          1205          5654
 9 SI_EK       35       1968       1968          1396          6523
10 SI_HA       41       1967       1967          2991         15143
11 TH           5       2000       2000            85           340

sample_size |>
  ungroup() |>
  summarise(sum_ind = sum(n_individuals), sum_n = sum(n_data_points))

ungroup: no grouping variables
summarise: now one row and 2 columns, ungrouped

# A tibble: 1 × 2
  sum_ind  sum_n
    <int>  <int>
1   37680 198323

write_csv(sample_size, paste0(home, "/output/sample_size.csv"))

# Compare how the means and quantiles differ depending on this filtering
IVBG_filter <- IVBG |>
  drop_na(k_se) |>
  #filter(k_se < quantile(k_se, probs = 0.95)) |> 
  filter(k_se < k)

drop_na: no rows removed
filter: removed 10,208 rows (5%), 188,115 rows remaining

# Summarize growth coefficients by cohort and area
VBG <- IVBG |>
  filter(k_se < k) |> # new!
  group_by(cohort, area) |>
  summarize(k_mean = mean(k, na.rm = T),
            k_median = quantile(k, prob = 0.5, na.rm = T),
            linf_median = quantile(linf, prob = 0.5, na.rm = T),
            k_lwr = quantile(k, prob = 0.05, na.rm = T),
            k_upr = quantile(k, prob = 0.95, na.rm = T)) |> 
  ungroup()

filter: removed 10,208 rows (5%), 188,115 rows remaining
group_by: 2 grouping variables (cohort, area)
summarize: now 343 rows and 7 columns, one group variable remaining (cohort)
ungroup: no grouping variables

VBG_filter <- IVBG_filter |>
  group_by(cohort, area) |>
  summarize(k_mean = mean(k, na.rm = T),
            k_median = quantile(k, prob = 0.5, na.rm = T),
            k_lwr = quantile(k, prob = 0.05, na.rm = T),
            k_upr = quantile(k, prob = 0.95, na.rm = T)) |> 
  ungroup()

group_by: 2 grouping variables (cohort, area)
summarize: now 343 rows and 6 columns, one group variable remaining (cohort)
ungroup: no grouping variables

ggplot() +
  geom_ribbon(data = VBG, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr,
                              fill = "All k's"), alpha = 0.4) +
  geom_ribbon(data = VBG_filter, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr,
                                     fill = "Filtered"), alpha = 0.4) +
  geom_line(data = VBG, aes(cohort, k_median, color = "All k's")) + 
  geom_line(data = VBG_filter, aes(cohort, k_median, color = "Filtered")) + 
  guides(fill = "none") +
  facet_wrap(~area)

ggplot() +
  geom_line(data = VBG, aes(cohort, k_median, color = "All k's")) + 
  geom_line(data = VBG_filter, aes(cohort, k_median, color = "Filtered")) + 
  guides(fill = "none") +
  facet_wrap(~area)

# No difference at all really. We should probably just discuss that with this model, achieving biologically reasonable parameter values and a good fit to data are sometimes two different things. In our case, we just want a representative value of the growth (as in time to reach average maximum size in the population) that accounts for the entire growth trajectory of an individual, for each area and cohort.

BT 1996 has a very high k, so I’ll inspect it in more detail

#|cache: false

bt96 <- IVBG |> 
  filter(area == "BT" & cohort %in% c(1995, 1996, 1997, 1998))

filter: removed 197,958 rows (>99%), 365 rows remaining

ggplot(bt96, aes(as.factor(cohort), k)) + 
  geom_jitter(width = 0.2, height = 0) + 
  geom_boxplot(fill = NA, width = 0.2) + 
  geom_hline(yintercept = mean(filter(bt96, cohort == 1996)$k), linetype = 2)

filter: removed 293 rows (80%), 72 rows remaining

# Check the fits in more detail
bt96b <- bt96 |> filter(cohort == 1996) |> as.data.frame()

filter: removed 293 rows (80%), 72 rows remaining

t <- bt96b |> 
  filter(k_se < k) |>
  filter(linf_se < linf) |> 
  filter(linf_lwr > 0)

filter: removed 16 rows (22%), 56 rows remaining

filter: no rows removed
filter: no rows removed

t

                        ID         k       k_se     linf    linf_se       k_lwr
1  2000_139_BT.female.4.NA 0.6998057 0.03282194 176.1500   3.041886  0.63547474
2  2000_139_BT.female.4.NA 0.6998057 0.03282194 176.1500   3.041886  0.63547474
3  2000_139_BT.female.4.NA 0.6998057 0.03282194 176.1500   3.041886  0.63547474
4  2000_139_BT.female.4.NA 0.6998057 0.03282194 176.1500   3.041886  0.63547474
5  2000_169_BT.female.4.NA 0.8309359 0.28906869 178.6083  20.099801  0.26436123
6  2000_169_BT.female.4.NA 0.8309359 0.28906869 178.6083  20.099801  0.26436123
7  2000_169_BT.female.4.NA 0.8309359 0.28906869 178.6083  20.099801  0.26436123
8  2000_169_BT.female.4.NA 0.8309359 0.28906869 178.6083  20.099801  0.26436123
9   2000_24_BT.female.4.NA 0.5728722 0.05217399 176.2271   6.819268  0.47061116
10  2000_24_BT.female.4.NA 0.5728722 0.05217399 176.2271   6.819268  0.47061116
11  2000_24_BT.female.4.NA 0.5728722 0.05217399 176.2271   6.819268  0.47061116
12  2000_24_BT.female.4.NA 0.5728722 0.05217399 176.2271   6.819268  0.47061116
13   2000_4_BT.female.4.NA 0.2457892 0.06377208 435.4096  75.698097  0.12079595
14   2000_4_BT.female.4.NA 0.2457892 0.06377208 435.4096  75.698097  0.12079595
15   2000_4_BT.female.4.NA 0.2457892 0.06377208 435.4096  75.698097  0.12079595
16   2000_4_BT.female.4.NA 0.2457892 0.06377208 435.4096  75.698097  0.12079595
17  2000_50_BT.female.4.NA 0.5000396 0.13927044 318.8768  41.323942  0.22706959
18  2000_50_BT.female.4.NA 0.5000396 0.13927044 318.8768  41.323942  0.22706959
19  2000_50_BT.female.4.NA 0.5000396 0.13927044 318.8768  41.323942  0.22706959
20  2000_50_BT.female.4.NA 0.5000396 0.13927044 318.8768  41.323942  0.22706959
21  2000_51_BT.female.4.NA 0.6280041 0.19388224 206.1642  25.353579  0.24799490
22  2000_51_BT.female.4.NA 0.6280041 0.19388224 206.1642  25.353579  0.24799490
23  2000_51_BT.female.4.NA 0.6280041 0.19388224 206.1642  25.353579  0.24799490
24  2000_51_BT.female.4.NA 0.6280041 0.19388224 206.1642  25.353579  0.24799490
25  2000_57_BT.female.4.NA 0.3710587 0.19779509 333.8080  98.764561 -0.01661973
26  2000_57_BT.female.4.NA 0.3710587 0.19779509 333.8080  98.764561 -0.01661973
27  2000_57_BT.female.4.NA 0.3710587 0.19779509 333.8080  98.764561 -0.01661973
28  2000_57_BT.female.4.NA 0.3710587 0.19779509 333.8080  98.764561 -0.01661973
29  2000_78_BT.female.4.NA 0.4577937 0.05308862 208.3834  11.887232  0.35374004
30  2000_78_BT.female.4.NA 0.4577937 0.05308862 208.3834  11.887232  0.35374004
31  2000_78_BT.female.4.NA 0.4577937 0.05308862 208.3834  11.887232  0.35374004
32  2000_78_BT.female.4.NA 0.4577937 0.05308862 208.3834  11.887232  0.35374004
33  2000_82_BT.female.4.NA 0.5529989 0.06526315 179.5668   9.224169  0.42508312
34  2000_82_BT.female.4.NA 0.5529989 0.06526315 179.5668   9.224169  0.42508312
35  2000_82_BT.female.4.NA 0.5529989 0.06526315 179.5668   9.224169  0.42508312
36  2000_82_BT.female.4.NA 0.5529989 0.06526315 179.5668   9.224169  0.42508312
37  2000_83_BT.female.4.NA 0.6063114 0.05190447 178.8868   6.254236  0.50457860
38  2000_83_BT.female.4.NA 0.6063114 0.05190447 178.8868   6.254236  0.50457860
39  2000_83_BT.female.4.NA 0.6063114 0.05190447 178.8868   6.254236  0.50457860
40  2000_83_BT.female.4.NA 0.6063114 0.05190447 178.8868   6.254236  0.50457860
41  2000_95_BT.female.4.NA 0.1261268 0.07828577 774.2764 390.036643 -0.02731335
42  2000_95_BT.female.4.NA 0.1261268 0.07828577 774.2764 390.036643 -0.02731335
43  2000_95_BT.female.4.NA 0.1261268 0.07828577 774.2764 390.036643 -0.02731335
44  2000_95_BT.female.4.NA 0.1261268 0.07828577 774.2764 390.036643 -0.02731335
45 2001_201_BT.female.5.NA 0.5572245 0.11202672 366.7864  31.931152  0.33765214
46 2001_201_BT.female.5.NA 0.5572245 0.11202672 366.7864  31.931152  0.33765214
47 2001_201_BT.female.5.NA 0.5572245 0.11202672 366.7864  31.931152  0.33765214
48 2001_201_BT.female.5.NA 0.5572245 0.11202672 366.7864  31.931152  0.33765214
49 2001_222_BT.female.5.NA 0.4204646 0.14090368 397.1261  68.884483  0.14429342
50 2001_222_BT.female.5.NA 0.4204646 0.14090368 397.1261  68.884483  0.14429342
51 2001_222_BT.female.5.NA 0.4204646 0.14090368 397.1261  68.884483  0.14429342
52 2001_222_BT.female.5.NA 0.4204646 0.14090368 397.1261  68.884483  0.14429342
53 2001_236_BT.female.5.NA 0.2093150 0.09910580 573.1898 192.597276  0.01506767
54 2001_236_BT.female.5.NA 0.2093150 0.09910580 573.1898 192.597276  0.01506767
55 2001_236_BT.female.5.NA 0.2093150 0.09910580 573.1898 192.597276  0.01506767
56 2001_236_BT.female.5.NA 0.2093150 0.09910580 573.1898 192.597276  0.01506767
       k_upr   linf_lwr  linf_upr row_id area_cohort_age area cohort
1  0.7641367 170.187889  182.1121  25583       BT 1996 1   BT   1996
2  0.7641367 170.187889  182.1121  25583       BT 1996 2   BT   1996
3  0.7641367 170.187889  182.1121  25583       BT 1996 3   BT   1996
4  0.7641367 170.187889  182.1121  25583       BT 1996 4   BT   1996
5  1.3975105 139.212721  218.0039  25667       BT 1996 1   BT   1996
6  1.3975105 139.212721  218.0039  25667       BT 1996 2   BT   1996
7  1.3975105 139.212721  218.0039  25667       BT 1996 3   BT   1996
8  1.3975105 139.212721  218.0039  25667       BT 1996 4   BT   1996
9  0.6751332 162.861337  189.5929  25905       BT 1996 1   BT   1996
10 0.6751332 162.861337  189.5929  25905       BT 1996 2   BT   1996
11 0.6751332 162.861337  189.5929  25905       BT 1996 3   BT   1996
12 0.6751332 162.861337  189.5929  25905       BT 1996 4   BT   1996
13 0.3707825 287.041290  583.7778  26101       BT 1996 1   BT   1996
14 0.3707825 287.041290  583.7778  26101       BT 1996 2   BT   1996
15 0.3707825 287.041290  583.7778  26101       BT 1996 3   BT   1996
16 0.3707825 287.041290  583.7778  26101       BT 1996 4   BT   1996
17 0.7730097 237.881859  399.8717  26103       BT 1996 1   BT   1996
18 0.7730097 237.881859  399.8717  26103       BT 1996 2   BT   1996
19 0.7730097 237.881859  399.8717  26103       BT 1996 3   BT   1996
20 0.7730097 237.881859  399.8717  26103       BT 1996 4   BT   1996
21 1.0080133 156.471211  255.8572  26104       BT 1996 1   BT   1996
22 1.0080133 156.471211  255.8572  26104       BT 1996 2   BT   1996
23 1.0080133 156.471211  255.8572  26104       BT 1996 3   BT   1996
24 1.0080133 156.471211  255.8572  26104       BT 1996 4   BT   1996
25 0.7587370 140.229501  527.3866  26119       BT 1996 1   BT   1996
26 0.7587370 140.229501  527.3866  26119       BT 1996 2   BT   1996
27 0.7587370 140.229501  527.3866  26119       BT 1996 3   BT   1996
28 0.7587370 140.229501  527.3866  26119       BT 1996 4   BT   1996
29 0.5618474 185.084375  231.6823  26157       BT 1996 1   BT   1996
30 0.5618474 185.084375  231.6823  26157       BT 1996 2   BT   1996
31 0.5618474 185.084375  231.6823  26157       BT 1996 3   BT   1996
32 0.5618474 185.084375  231.6823  26157       BT 1996 4   BT   1996
33 0.6809147 161.487402  197.6461  26163       BT 1996 1   BT   1996
34 0.6809147 161.487402  197.6461  26163       BT 1996 2   BT   1996
35 0.6809147 161.487402  197.6461  26163       BT 1996 3   BT   1996
36 0.6809147 161.487402  197.6461  26163       BT 1996 4   BT   1996
37 0.7080441 166.628541  191.1451  26166       BT 1996 1   BT   1996
38 0.7080441 166.628541  191.1451  26166       BT 1996 2   BT   1996
39 0.7080441 166.628541  191.1451  26166       BT 1996 3   BT   1996
40 0.7080441 166.628541  191.1451  26166       BT 1996 4   BT   1996
41 0.2795669   9.804558 1538.7482  26194       BT 1996 1   BT   1996
42 0.2795669   9.804558 1538.7482  26194       BT 1996 2   BT   1996
43 0.2795669   9.804558 1538.7482  26194       BT 1996 3   BT   1996
44 0.2795669   9.804558 1538.7482  26194       BT 1996 4   BT   1996
45 0.7767969 304.201370  429.3715  26927       BT 1996 1   BT   1996
46 0.7767969 304.201370  429.3715  26927       BT 1996 2   BT   1996
47 0.7767969 304.201370  429.3715  26927       BT 1996 3   BT   1996
48 0.7767969 304.201370  429.3715  26927       BT 1996 4   BT   1996
49 0.6966359 262.112502  532.1397  27021       BT 1996 1   BT   1996
50 0.6966359 262.112502  532.1397  27021       BT 1996 2   BT   1996
51 0.6966359 262.112502  532.1397  27021       BT 1996 3   BT   1996
52 0.6966359 262.112502  532.1397  27021       BT 1996 4   BT   1996
53 0.4035624 195.699170  950.6805  27079       BT 1996 1   BT   1996
54 0.4035624 195.699170  950.6805  27079       BT 1996 2   BT   1996
55 0.4035624 195.699170  950.6805  27079       BT 1996 3   BT   1996
56 0.4035624 195.699170  950.6805  27079       BT 1996 4   BT   1996

bt96 |> 
  filter(k_se < k) |>
  filter(linf_se < linf) |> 
  filter(linf_lwr > 0) |> 
  ggplot(aes(as.factor(cohort), k)) + 
  geom_jitter(width = 0.2, height = 0) + 
  geom_boxplot(fill = NA, width = 0.2) + 
  geom_hline(yintercept = mean(filter(bt96, cohort == 1996)$k), linetype = 2)

filter: removed 45 rows (12%), 320 rows remaining
filter: no rows removed

filter: removed 33 rows (10%), 287 rows remaining

filter: removed 293 rows (80%), 72 rows remaining

# Plot these individuals in t
t <- t |> 
  dplyr::select(-cohort) |> 
  left_join(d, by = "ID") |> 
  mutate(length_mm_pred = linf*(1-exp(-k*age_bc)))

left_join: added 15 columns (area_cohort_age.x, area.x, length_mm, age_bc, age_catch, …)

           > rows only in x         0

           > rows only in y  (295,745)

           > matched rows         224    (includes duplicates)

           >                 =========

           > rows total           224

mutate: new variable 'length_mm_pred' (double) with 56 unique values and 0% NA

ggplot(t, aes(as.factor(cohort), k)) + 
  geom_jitter(width = 0.2, height = 0) + 
  geom_boxplot(fill = NA, width = 0.2) + 
  geom_hline(yintercept = mean(t$k), linetype = 2)

ggplot(t, aes(age_bc, length_mm)) +
  geom_jitter(width = 0.3) +
  geom_line(aes(age_bc, length_mm_pred)) +
  facet_wrap(~ID)

# These fits look totally reasonable

Add GAM-predicted temperature to growth models

pred_temp <- read_csv(paste0(home, "/output/gam_predicted_temps.csv")) |> 
  rename(cohort = year)

Rows: 395 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): area, model
dbl (2): year, temp

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed one variable (cohort)

VBG <- VBG |>
  left_join(pred_temp, by = c("area", "cohort"))

left_join: added 2 columns (temp, model)
           > rows only in x     0
           > rows only in y  ( 52)
           > matched rows     343
           >                 =====
           > rows total       343

# Save data for map-plot
cohort_sample_size <- IVBG |>
  group_by(area, cohort) |> 
  summarise(n = n()) # individuals, not samples!

group_by: 2 grouping variables (area, cohort)
summarise: now 343 rows and 3 columns, one group variable remaining (area)

VBG <- left_join(VBG, cohort_sample_size, by = c("cohort", "area"))

left_join: added one column (n)
           > rows only in x     0
           > rows only in y  (  0)
           > matched rows     343
           >                 =====
           > rows total       343

write_csv(VBG, paste0(home, "/output/vbg.csv"))

# Calculate the plotting order (also used for map plot)
order <- VBG |>
  ungroup() |>
  group_by(area) |>
  summarise(min_temp = min(temp)) |>
  arrange(desc(min_temp))

ungroup: no grouping variables
group_by: one grouping variable (area)
summarise: now 11 rows and 2 columns, ungrouped

order

# A tibble: 11 × 2
   area  min_temp
   <chr>    <dbl>
 1 SI_HA     7.83
 2 TH        6.85
 3 BS        6.22
 4 BT        5.87
 5 FM        5.87
 6 SI_EK     5.49
 7 MU        5.16
 8 FB        5.04
 9 HO        3.99
10 JM        3.37
11 RA        3.06

nareas <- length(unique(order$area)) + 2 # to skip the brightest colors
colors <- colorRampPalette(brewer.pal(name = "RdYlBu", n = 10))(nareas)[-c(7,8)]

Plot VBGE fits

# Sample 50 IDs per area and plot their data and VBGE fits
set.seed(4)
ids <- IVBG |> distinct(ID, .keep_all = TRUE) |> group_by(area) |> slice_sample(n = 30)

distinct: removed 160,643 rows (81%), 37,680 rows remaining

group_by: one grouping variable (area)

slice_sample (grouped): removed 37,350 rows (99%), 330 rows remaining

#ids |> ungroup() |> group_by(area) |> summarise(n = length(unique(ID))) |> arrange(n)

IVBG2 <- IVBG |>
  filter(ID %in% ids$ID) |> 
  distinct(ID, .keep_all = TRUE) |> 
  select(ID, k, linf)

filter: removed 196,691 rows (99%), 1,632 rows remaining

distinct: removed 1,302 rows (80%), 330 rows remaining

select: dropped 10 variables (k_se, linf_se, k_lwr, k_upr, linf_lwr, …)

d2 <- d |>
  ungroup() |>
  filter(ID %in% ids$ID) |>
  left_join(IVBG2, by = "ID") |>
  mutate(length_mm_pred = linf*(1-exp(-k*age_bc)))

ungroup: no grouping variables

filter: removed 294,169 rows (99%), 1,632 rows remaining

left_join: added 2 columns (k, linf)

           > rows only in x       0

           > rows only in y  (    0)

           > matched rows     1,632

           >                 =======

           > rows total       1,632

mutate: new variable 'length_mm_pred' (double) with 1,632 unique values and 0% NA

ggplot(d2, aes(age_bc, length_mm, group = ID, color = ID)) +
  geom_jitter(width = 0.3, height = 0, alpha = 0.5, size = 0.4) +
  geom_line(aes(age_bc, length_mm_pred, group = ID, color = ID),
            inherit.aes = FALSE, alpha = 0.8, linewidth = 0.3) + 
  guides(color = "none") +
  scale_color_viridis(discrete = TRUE, name = "Area", option = "cividis") +
  labs(x = "Age (yrs)", y = "Length (mm)") +
  facet_wrap(~factor(area, levels = (arrange(area_attr, desc(lat)))$area), ncol = 3)

ggsave(paste0(home, "/figures/vb_fits.pdf" ), width = 17, height = 17, units = "cm")

# rugs <- IVBG |> 
#   group_by(area) |> 
#   summarise(median_k = median(k),
#             median_linf = median(linf))

# k <- IVBG |> 
#   ggplot(aes(k, color = factor(area, order$area))) + 
#   geom_density(alpha = 0.4, fill = NA, adjust = 1.5) +  
#   scale_color_manual(values = colors, name = "Area") +
#   coord_cartesian(expand = 0) + 
#   labs(x = expression(italic(k))) +
#   geom_rug(data = rugs, aes(median_k)) +
#   theme(legend.position = c(0.9, 0.5))

# violons instead?
k2 <- IVBG |> 
  ggplot(aes(factor(area, order$area), k, fill = factor(area, order$area))) + 
  geom_violin(alpha = 0.6, color = NA) +  
  geom_boxplot(outlier.color = NA, width = 0.2, alpha = 0.8, fill = NA, size = 0.4) +
  # geom_jitter(data = IVBG |> group_by(area) |> slice_sample(n = 100),
  #             aes(factor(area, order$area), k), height = 0, width = 0.2, alpha = 0.1,
  #             color = "grey40") +
  scale_fill_manual(values = colors, name = "Area") +
  guides(fill = "none", color = "none") +
  labs(x = "", y = expression(italic(k)))

k2

# linf <- IVBG |> 
#   ggplot(aes(linf, color = factor(area, order$area))) + 
#   geom_density(alpha = 0.4, fill = NA, adjust = 1.5) +  
#   scale_color_manual(values = colors, name = "Area") +
#   coord_cartesian(expand = 0, xlim = c(0, 2000)) + 
#   labs(x = expression(paste(italic(L[infinity]), " [cm]"))) +
#   geom_rug(data = rugs, aes(median_linf)) +
#   guides(color = "none")

linf2 <- IVBG |> 
  # group_by(area) |> 
  # mutate(upr = quantile(linf, probs = 0.95)) |> 
  # ungroup() |> 
  # filter(linf < upr) |> 
  ggplot(aes(factor(area, order$area), linf, fill = factor(area, order$area))) + 
  geom_violin(alpha = 0.6, color = NA) +  
  geom_boxplot(outlier.color = NA, width = 0.1, alpha = 0.8, fill = NA, size = 0.4) +
  coord_cartesian(ylim = c(0, 2300)) +
  # geom_jitter(data = IVBG |> group_by(area) |> slice_sample(n = 100),
  #             aes(factor(area, order$area), linf), height = 0, width = 0.2, alpha = 0.1,
  #             color = "grey40") +
  scale_fill_manual(values = colors, name = "Area") +
  guides(fill = "none", color = "none") +
  labs(x = "", y = expression(paste(italic(L[infinity]), " [cm]")))

linf2

#k / linf
k2 / linf2

ggsave(paste0(home, "/figures/vb_pars.pdf" ), width = 17, height = 17, units = "cm")

Fit Sharpe-Schoolfield model to K

By area

model <- 'sharpeschoolhigh_1981'

# Get starting values on full dataset for Sharpe-Schoolfield model
dat <- VBG |>
  select(k_median, temp, area) |>
  rename(rate = k_median)

select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)

rename: renamed one variable (rate)

lower <- get_lower_lims(dat$temp, dat$rate, model_name = model)
upper <- get_upper_lims(dat$temp, dat$rate, model_name = model)
start <- get_start_vals(dat$temp, dat$rate, model_name = model)
  
# Sharpe-Schoolfield model fit to data for each area
preds <- NULL
pred <- NULL
pars <- list()

for(a in unique(dat$area)) {
  
  # Get data
  dd <- dat |> filter(area == a)
  
  # Fit model
  fit <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dd,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )
  
  # Make predictions on new data
  new_data <- data.frame(temp = seq(min(dd$temp), max(dd$temp), length.out = 100))
  
  pred <- augment(fit, newdata = new_data) |> mutate(area = a)
  
  # Add to general data frame
  preds <- data.frame(rbind(preds, pred))
  
  # Add parameter estimates
  pars[[a]] <- summary(fit)$coefficients |>
    as.data.frame() |>
    rownames_to_column(var = "parameter") |>
    mutate(area = a)
  
}

filter: removed 280 rows (82%), 63 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 305 rows (89%), 38 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 302 rows (88%), 41 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 308 rows (90%), 35 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 295 rows (86%), 48 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 316 rows (92%), 27 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 324 rows (94%), 19 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 309 rows (90%), 34 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 325 rows (95%), 18 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 328 rows (96%), 15 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

filter: removed 338 rows (99%), 5 rows remaining

mutate: new variable 'area' (character) with one unique value and 0% NA
mutate: new variable 'area' (character) with one unique value and 0% NA

# For comparison with the bayesian model below
nls_pars <- bind_rows(pars)

nls_pars <- nls_pars |> 
  mutate(upr = Estimate + 1.96 *`Std. Error`,
         lwr = Estimate - 1.96 *`Std. Error`) |> 
  dplyr::select(parameter, Estimate, upr, lwr, area) |> 
  mutate(source = "nls")

mutate: new variable 'upr' (double) with 44 unique values and 0% NA

        new variable 'lwr' (double) with 44 unique values and 0% NA

mutate: new variable 'source' (character) with one unique value and 0% NA

All areas pooled

# One sharpe schoolfield fitted to all data
fit_all <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dat,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )

# Make predictions on new data
new_data_all <- data.frame(temp = seq(min(dat$temp), max(dat$temp), length.out = 100))

pred_all <- augment(fit_all, newdata = new_data_all) |> 
  mutate(area = "all")

mutate: new variable 'area' (character) with one unique value and 0% NA

# Add t_opt (not correct equation I think!) from Padfield 2021 ISME
kb <- 8.62e-05
# data.frame(par = names(coef(fit_all)), est = coef(fit_all)) |> 
#   pivot_wider(names_from = par, values_from = est) |> 
#   summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))
# 
# get_topt

# This rTPC function just finds the temperature where the function is maximized,
# do that also for brms fits
get_topt(fit_all)

[1] 10.00058

Plot Sharpe Schoolfield fits

# Plot growth coefficients by year and area against mean SST

p1 <- preds |>
  ggplot(aes(temp, .fitted, color = factor(area, order$area))) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 1, alpha = 0.8) +
  geom_line(linewidth = 1) +
  geom_line(data = pred_all, aes(temp, .fitted), linewidth = 1, inherit.aes = FALSE, linetype = 2) +
  scale_color_manual(values = colors, name = "Area") +
  guides(color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  scale_x_continuous(breaks = seq(-5, 20, 1)) +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)") +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal")

p1 + facet_wrap(~factor(area, levels = (arrange(area_attr, desc(lat)))$area))

p1

ggsave(paste0(home, "/figures/sharpe_school.pdf" ), width = 17, height = 17, units = "cm")

Can we fit a single Sharpe Schoolfield with area-specific parameters with brms?

# Again, here are the data we are fitting:
ggplot(dat, aes(temp, rate, color = factor(area, order$area))) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6)  

# Which family to use? Well we can try with a Gaussian or a student-t first... 
# After that, maybe gamma or lognormal, if we predict negative K's
hist(dat$rate)

# Here's the equation:
sharpeschoolhigh_1981

function (temp, r_tref, e, eh, th, tref) 
{
    tref <- 273.15 + tref
    k <- 8.62e-05
    boltzmann.term <- r_tref * exp(e/k * (1/tref - 1/(temp + 
        273.15)))
    inactivation.term <- 1/(1 + exp(eh/k * (1/(th + 273.15) - 
        1/(temp + 273.15))))
    return(boltzmann.term * inactivation.term)
}
<bytecode: 0x7f9d95f8aa30>
<environment: namespace:rTPC>

# Add in fixed parameters
dat$bk <- 8.62e-05
dat$tref <- 8 + 273.15

# We can use the nls() parameters as starting values
nlspars <- summary(fit_all)$coefficients[, 1]
summary(fit_all)

Formula: rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, 
    tref = 8)

Parameters:
       Estimate Std. Error t value Pr(>|t|)  
r_tref   0.5265     0.5818   0.905   0.3661  
e        1.3179     1.1940   1.104   0.2705  
eh       1.8364     0.7329   2.506   0.0127 *
th       6.5328     6.1852   1.056   0.2916  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.05265 on 339 degrees of freedom

Number of iterations to convergence: 33 
Achieved convergence tolerance: 1.49e-08

# Parameters:
#        Estimate Std. Error t value Pr(>|t|)  
# r_tref   0.5265     0.5818   0.905   0.3661  
# e        1.3179     1.1940   1.104   0.2705  
# eh       1.8364     0.7329   2.506   0.0127 *
# th       6.5328     6.1852   1.056   0.2916  

# As a first guess, I use the nls estimate as the mean, and standard deviation something close to it. We can also bound som fo them by 0, since negative values shouldn't be possible
# (better visualization including bounds further down)
n=10000
hist(rnorm(0.5, 2.5, n=n), main = "rtref")

hist(rnorm(1, 2.5, n=n), main = "e") 

hist(rnorm(2, 2.5, n=n), main = "eh")

hist(rnorm(6, 5, n=n), main = "hh")

fit_prior <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1, 
     nl = TRUE),
  data = dat,
  cores = 2,
  chains = 2,
  iter = 1500,
  sample_prior = "only",
  seed = 9,
  prior = c(
    prior(normal(0.5, 2.5), nlpar = "rtref", lb = 0),
    prior(normal(1, 2.5), nlpar = "e", lb = 0),
    prior(normal(2, 2.5), nlpar = "eh", lb = 0),
    prior(normal(6, 5), nlpar = "th", lb = 0)
  )
)

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

# Global prior predictive check in relation to data. Doesn't loo too informative...
dat |>
  data_grid(temp = seq_range(temp, n = 51)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fit_prior) |> 
  mutate(rtref = nlspars[1], e = nlspars[2], eh = nlspars[3], th = nlspars[4]) |> 
  mutate(man_pred = (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))) |> 
  ggplot(aes(temp, y = .epred)) +
  stat_lineribbon(.width = c(0.95), alpha = 0.3, color = "black", fill = "black") + 
  geom_line(aes(y = man_pred), color = "tomato3", linetype = 2) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Area") +
  coord_cartesian(ylim = c(0, 10)) + # Extends waaaay higher... 
  NULL

ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

mutate (grouped): new variable 'rtref' (double) with one unique value and 0% NA

                  new variable 'e' (double) with one unique value and 0% NA

                  new variable 'eh' (double) with one unique value and 0% NA

                  new variable 'th' (double) with one unique value and 0% NA

mutate (grouped): new variable 'man_pred' (double) with 51 unique values and 0% NA

# Now make sure it converges with real data but no random effects
fit_global <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1,
     nl = TRUE),
  data = dat,
  cores = 4,
  chains = 4,
  iter = 4000, 
  seed = 9,
  sample_prior = "yes",
  prior = c(
    prior(normal(0.5, 2.5), nlpar = "rtref", lb = 0),
    prior(normal(1, 2.5), nlpar = "e", lb = 0),
    prior(normal(2, 2.5), nlpar = "eh", lb = 0),
    prior(normal(6, 5), nlpar = "th", lb = 0)
  ),
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

plot(fit_global)

pp_check(fit_global)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

plot(conditional_effects(fit_global, effect = "temp"), points = TRUE)

dat |>
  data_grid(temp = seq_range(temp, n = 51)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fit_global) |> 
  mutate(rtref = nlspars[1], e = nlspars[2], eh = nlspars[3], th = nlspars[4]) |> # from summary(fit_all)
  mutate(man_pred = (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15))))) |> 
  ggplot(aes(temp, y = .epred)) +
  stat_lineribbon(.width = c(0.95), alpha = 0.3, color = "black", fill = "black") + 
  geom_line(aes(y = man_pred), color = "tomato3", linetype = 2) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Area")

ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

mutate (grouped): new variable 'rtref' (double) with one unique value and 0% NA

                  new variable 'e' (double) with one unique value and 0% NA

                  new variable 'eh' (double) with one unique value and 0% NA

                  new variable 'th' (double) with one unique value and 0% NA

mutate (grouped): new variable 'man_pred' (double) with 51 unique values and 0% NA

# Plot prior vs posterior
post <- fit_global |>
  as_draws_df() |>
  dplyr::select(b_rtref_Intercept, b_e_Intercept, b_eh_Intercept, b_th_Intercept) |> 
  rename(rtref = b_rtref_Intercept, 
         e = b_e_Intercept,
         eh = b_eh_Intercept,
         th = b_th_Intercept) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Posterior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (rtref, e, eh, th)

pivot_longer: reorganized (rtref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

prior <- fit_global |>
  as_draws_df() |>
  dplyr::select(prior_b_rtref, prior_b_e, prior_b_eh, prior_b_th) |> 
    rename(rtref = prior_b_rtref, 
         e = prior_b_e,
         eh = prior_b_eh,
         th = prior_b_th) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Prior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (rtref, e, eh, th)

pivot_longer: reorganized (rtref, e, eh, th) into (parameter, value) [was 8000x4, now 32000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

dist <- bind_rows(prior, post)

ggplot(dist, aes(value, fill = type)) +
  geom_density(color = NA, alpha = 0.5) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~parameter, scales = "free") + 
  theme(legend.position = c(0.35, 0.9)) +
  labs(fill = "")

Ok, seems to work with priors and everything. They are roughly as broad as can be and still have a fit. Now fit the full model, with random area effects, more iterations and chains!

# Still good enough, might change some priors. Now look at random area effects!
# Gaussian model
fitb <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1 + (1|area),
     nl = TRUE),
  data = dat,
  cores = 4,
  chains = 4,
  iter = 5000,
  sample_prior = "yes",
  save_pars = save_pars(all = TRUE),
  seed = 9,
  prior = c(
    prior(normal(0.5, 2.5), nlpar = "rtref", lb = 0),
    prior(normal(1, 2.5), nlpar = "e", lb = 0),
    prior(normal(2, 2.5), nlpar = "eh", lb = 0),
    prior(normal(6, 5), nlpar = "th", lb = 0)
  ),
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

# Student-t model
fitbs <- brm(
  bf(rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15)))) / (1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))),
     rtref + e + eh + th ~ 1 + (1|area),
     nl = TRUE),
  data = dat,
  family = student,
  cores = 4,
  chains = 4,
  iter = 5000,
  sample_prior = "yes",
  save_pars = save_pars(all = TRUE),
  seed = 9,
  prior = c(
    prior(normal(0.5, 2.5), nlpar = "rtref", lb = 0),
    prior(normal(1, 2.5), nlpar = "e", lb = 0),
    prior(normal(2, 2.5), nlpar = "eh", lb = 0),
    prior(normal(6, 5), nlpar = "th", lb = 0)
  ),
  control = list(adapt_delta = 0.99, max_treedepth = 12)
)

Compiling Stan program...
Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.2/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DBOOST_NO_AUTO_PTR  -include '/Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/usr/local/include   -fPIC  -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:88:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
namespace Eigen {
^
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
namespace Eigen {
               ^
               ;
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
In file included from /Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Dense:1:
/Library/Frameworks/R.framework/Versions/4.2/Resources/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#include <complex>
         ^~~~~~~~~
3 errors generated.
make: *** [foo.o] Error 1

Start sampling

Compare fitted models based on ELPD. Preferred model in the first row!

fitb_loo <- loo(fitb, moment_match = TRUE)

Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Warning: Found 1 observations with a pareto_k > 0.7 in model 'fitb'. It is
recommended to set 'reloo = TRUE' in order to calculate the ELPD without the
assumption that these observations are negligible. This will refit the model 1
times to compute the ELPDs for the problematic observations directly.

fitbs_loo <- loo(fitbs, moment_match = TRUE)

Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Warning: Found 2 observations with a pareto_k > 0.7 in model 'fitbs'. It is
recommended to set 'reloo = TRUE' in order to calculate the ELPD without the
assumption that these observations are negligible. This will refit the model 2
times to compute the ELPDs for the problematic observations directly.

loo_compare(fitb_loo, fitbs_loo)

      elpd_diff se_diff
fitbs   0.0       0.0  
fitb  -10.8      10.8  

The student model is preferred. Inspect it!

# Check fit
summary(fitbs)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))) 
         rtref ~ 1 + (1 | area)
         e ~ 1 + (1 | area)
         eh ~ 1 + (1 | area)
         th ~ 1 + (1 | area)
   Data: dat (Number of observations: 343) 
  Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~area (Number of levels: 11) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(rtref_Intercept)     0.14      0.16     0.00     0.58 1.00     2574     4177
sd(e_Intercept)         0.21      0.13     0.01     0.52 1.00     2850     3740
sd(eh_Intercept)        0.16      0.18     0.01     0.64 1.00     2425     2599
sd(th_Intercept)        0.48      0.57     0.02     2.02 1.00     2514     3850

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept     1.55      1.24     0.30     4.78 1.00     1878     2279
e_Intercept         1.89      0.93     0.36     3.71 1.00     1931     2368
eh_Intercept        2.21      0.83     0.67     3.89 1.00     2309     2601
th_Intercept        4.43      3.05     0.56    12.77 1.00     1863     2153

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.04      0.00     0.03     0.04 1.00    12112     8196
nu       10.37      4.30     5.10    21.61 1.00    11772     7771

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(fitbs)

pp_check(fitbs)

Using 10 posterior draws for ppc type 'dens_overlay' by default.

ggsave(paste0(home, "/figures/supp/sharpe_brms_ppcheck.pdf" ), width = 11, height = 11, units = "cm")

# Plot prior vs posterior
post <- fitbs |>
  as_draws_df() |>
  dplyr::select(b_rtref_Intercept, b_e_Intercept, b_eh_Intercept, b_th_Intercept) |> 
  rename(rtref = b_rtref_Intercept, 
         e = b_e_Intercept,
         eh = b_eh_Intercept,
         th = b_th_Intercept) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Posterior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (rtref, e, eh, th)

pivot_longer: reorganized (rtref, e, eh, th) into (parameter, value) [was 10000x4, now 40000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

prior <- fitbs |>
  as_draws_df() |>
  dplyr::select(prior_b_rtref, prior_b_e, prior_b_eh, prior_b_th) |> 
    rename(rtref = prior_b_rtref, 
         e = prior_b_e,
         eh = prior_b_eh,
         th = prior_b_th) |> 
  pivot_longer(everything(), names_to = "parameter") |> 
  mutate(type = "Prior")

Warning: Dropping 'draws_df' class as required metadata was removed.

rename: renamed 4 variables (rtref, e, eh, th)

pivot_longer: reorganized (rtref, e, eh, th) into (parameter, value) [was 10000x4, now 40000x2]

mutate: new variable 'type' (character) with one unique value and 0% NA

dist <- bind_rows(prior, post)

ggplot(dist, aes(value, fill = type)) +
  geom_density(color = NA, alpha = 0.5) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~parameter, scales = "free") + 
  theme(legend.position = c(0.35, 0.9)) +
  labs(fill = "")

ggsave(paste0(home, "/figures/supp/sharpe_brms_prior.pdf" ), width = 17, height = 17, units = "cm")

# Calculate T_opt
ts <- dat |> 
  data_grid(temp = seq_range(temp, n = 100)) |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs, re_formula = NA, ndraws = 10000) 

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ggplot(ts) + 
  geom_line(aes(temp, y = .epred, group = .draw), alpha = .1)

# Compute quantiles across spaghetties
t_opt <- ts |> 
  group_by(.draw) |> 
  filter(.epred == max(.epred)) |> 
  ungroup()

group_by: one grouping variable (.draw)

filter (grouped): removed 990,000 rows (99%), 10,000 rows remaining

ungroup: no grouping variables

ggplot(t_opt, aes(temp)) + 
  geom_histogram()

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

quantile(t_opt$temp, probs = c(0.025, 0.5, 0.975))

     2.5%       50%     97.5% 
 6.960337  9.468607 16.854070 

# Make the main plot (conditional effect of temperature, with and without random effects)
# Predictions without random effects
glob <- dat |>
  data_grid(temp = seq_range(temp, n = 100)) |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs, re_formula = NA) |> 
  ungroup()

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

dat |>
  group_by(area) |>
  data_grid(temp = seq_range(temp, n = 100)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs) |> 
  ungroup() |> 
  ggplot(aes(temp, y = .epred, color = factor(area, order$area))) +
  stat_lineribbon(data = glob, aes(temp, .epred), .width = c(0.9), inherit.aes = FALSE,
                  fill = "black", color = "black", alpha = 0.1) +
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 1, alpha = 0.8) +
  stat_lineribbon(.width = c(0)) +
  stat_lineribbon(data = glob, aes(temp, .epred), .width = c(0), inherit.aes = FALSE,
                  color = "black", alpha = 0.9, linetype = 2) +
  # trimming the y-axis a bit... highest k is 0.53, second highest is 0.37
  coord_cartesian(ylim = c(min(dat$rate) - 0.01, sort(dat$rate)[length(dat$rate) - 1] + 0.01)) +
  scale_color_manual(values = colors, name = "Area") +
  guides(fill = "none",
         color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)")

group_by: one grouping variable (area)
ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

ungroup: no grouping variables

ggsave(paste0(home, "/figures/sharpe_school_bayes.pdf" ), width = 17, height = 17, units = "cm")

Now plot comparisons with the nls multstart fits

# Compare with area-specific sharpe scool!
area_pred_brms <- dat |>
  group_by(area) |> 
  data_grid(temp = seq_range(temp, n = 100)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15) |> 
  add_epred_draws(fitbs)

group_by: one grouping variable (area)

ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

# Plot area specific in comparison with nls multstart
p1 <- ggplot(area_pred_brms, aes(temp, y = .epred, color = factor(area, order$area),
                           fill = factor(area, order$area))) +
  stat_lineribbon(.width = c(0.95), alpha = 0.4, color = NA) +
  stat_lineribbon(.width = c(0)) + 
  geom_point(data = dat, aes(temp, rate, color = factor(area, order$area)), size = 0.6) +
  scale_color_manual(values = colors, name = "Area") +
  scale_fill_manual(values = colors, name = "Area") +
  scale_linetype_manual(values = 2) +
  facet_wrap(~factor(area, rev(order$area))) +
  guides(fill = "none",
         color = guide_legend(nrow = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1))) +
  theme(axis.title.y = ggtext::element_markdown(),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)",
       linetype = "")

p1

ggsave(paste0(home, "/figures/supp/sharpe_school_bayes_ci_facet.pdf" ), width = 17, height = 17, units = "cm")

# Now compare to nls fits
p1 + geom_line(data = preds, aes(temp, .fitted, linetype = "nls multstart"), color = "gray10")

# And for the population-level prediction
nd_brms <- dat |>
  data_grid(temp = seq_range(temp, n = 100)) |> 
  ungroup() |> 
  mutate(bk = 8.62e-05,
         tref = 8 + 273.15)

ungroup: no grouping variables

mutate: new variable 'bk' (double) with one unique value and 0% NA

        new variable 'tref' (double) with one unique value and 0% NA

pred_brms <- predict(fitbs, newdata = nd_brms, re_formula = NA) |> as.data.frame()
nd_brms$pred_pop <- pred_brms$Estimate

# (why not also comparing the global brms model...)
pred_brms_global <- predict(fit_global, newdata = nd_brms) |> as.data.frame()
nd_brms$pred_glob <- pred_brms_global$Estimate

ggplot(nd_brms, aes(temp, y = pred_pop, color = "mixed brms model")) +
  geom_point(data = dat, aes(temp, rate), size = 0.6, color = "grey30") +
  geom_line() + 
  geom_line(data = pred_all, aes(temp, .fitted, color = "nls"), linewidth = 1, inherit.aes = FALSE) +
  geom_line(data = nd_brms, aes(temp, pred_glob, color = "global brms model"), linewidth = 1, inherit.aes = FALSE) +
  guides(fill = "none",
         color = guide_legend(ncol = 1, reverse = TRUE, title.position = "top", title.hjust = 0.5,
                              override.aes = list(size = 1)),
         linetype = guide_legend(keywidth = unit(1, "cm"))) +
  theme(axis.title.y = ggtext::element_markdown()) +
  labs(x = "Temperature (°C)",
       y = "von Bertalanffy growth coefficient (*k*)",
       color = "")

And a comparison of parameter estimates…

get_variables(fitbs)

  [1] "b_rtref_Intercept"              "b_e_Intercept"                 
  [3] "b_eh_Intercept"                 "b_th_Intercept"                
  [5] "sd_area__rtref_Intercept"       "sd_area__e_Intercept"          
  [7] "sd_area__eh_Intercept"          "sd_area__th_Intercept"         
  [9] "sigma"                          "nu"                            
 [11] "r_area__rtref[BS,Intercept]"    "r_area__rtref[BT,Intercept]"   
 [13] "r_area__rtref[FB,Intercept]"    "r_area__rtref[FM,Intercept]"   
 [15] "r_area__rtref[HO,Intercept]"    "r_area__rtref[JM,Intercept]"   
 [17] "r_area__rtref[MU,Intercept]"    "r_area__rtref[RA,Intercept]"   
 [19] "r_area__rtref[SI_EK,Intercept]" "r_area__rtref[SI_HA,Intercept]"
 [21] "r_area__rtref[TH,Intercept]"    "r_area__e[BS,Intercept]"       
 [23] "r_area__e[BT,Intercept]"        "r_area__e[FB,Intercept]"       
 [25] "r_area__e[FM,Intercept]"        "r_area__e[HO,Intercept]"       
 [27] "r_area__e[JM,Intercept]"        "r_area__e[MU,Intercept]"       
 [29] "r_area__e[RA,Intercept]"        "r_area__e[SI_EK,Intercept]"    
 [31] "r_area__e[SI_HA,Intercept]"     "r_area__e[TH,Intercept]"       
 [33] "r_area__eh[BS,Intercept]"       "r_area__eh[BT,Intercept]"      
 [35] "r_area__eh[FB,Intercept]"       "r_area__eh[FM,Intercept]"      
 [37] "r_area__eh[HO,Intercept]"       "r_area__eh[JM,Intercept]"      
 [39] "r_area__eh[MU,Intercept]"       "r_area__eh[RA,Intercept]"      
 [41] "r_area__eh[SI_EK,Intercept]"    "r_area__eh[SI_HA,Intercept]"   
 [43] "r_area__eh[TH,Intercept]"       "r_area__th[BS,Intercept]"      
 [45] "r_area__th[BT,Intercept]"       "r_area__th[FB,Intercept]"      
 [47] "r_area__th[FM,Intercept]"       "r_area__th[HO,Intercept]"      
 [49] "r_area__th[JM,Intercept]"       "r_area__th[MU,Intercept]"      
 [51] "r_area__th[RA,Intercept]"       "r_area__th[SI_EK,Intercept]"   
 [53] "r_area__th[SI_HA,Intercept]"    "r_area__th[TH,Intercept]"      
 [55] "prior_b_rtref"                  "prior_b_e"                     
 [57] "prior_b_eh"                     "prior_b_th"                    
 [59] "prior_sigma"                    "prior_nu"                      
 [61] "prior_sd_area"                  "prior_sd_area__1"              
 [63] "prior_sd_area__2"               "prior_sd_area__3"              
 [65] "lprior"                         "lp__"                          
 [67] "z_1[1,1]"                       "z_1[1,2]"                      
 [69] "z_1[1,3]"                       "z_1[1,4]"                      
 [71] "z_1[1,5]"                       "z_1[1,6]"                      
 [73] "z_1[1,7]"                       "z_1[1,8]"                      
 [75] "z_1[1,9]"                       "z_1[1,10]"                     
 [77] "z_1[1,11]"                      "z_2[1,1]"                      
 [79] "z_2[1,2]"                       "z_2[1,3]"                      
 [81] "z_2[1,4]"                       "z_2[1,5]"                      
 [83] "z_2[1,6]"                       "z_2[1,7]"                      
 [85] "z_2[1,8]"                       "z_2[1,9]"                      
 [87] "z_2[1,10]"                      "z_2[1,11]"                     
 [89] "z_3[1,1]"                       "z_3[1,2]"                      
 [91] "z_3[1,3]"                       "z_3[1,4]"                      
 [93] "z_3[1,5]"                       "z_3[1,6]"                      
 [95] "z_3[1,7]"                       "z_3[1,8]"                      
 [97] "z_3[1,9]"                       "z_3[1,10]"                     
 [99] "z_3[1,11]"                      "z_4[1,1]"                      
[101] "z_4[1,2]"                       "z_4[1,3]"                      
[103] "z_4[1,4]"                       "z_4[1,5]"                      
[105] "z_4[1,6]"                       "z_4[1,7]"                      
[107] "z_4[1,8]"                       "z_4[1,9]"                      
[109] "z_4[1,10]"                      "z_4[1,11]"                     
[111] "accept_stat__"                  "stepsize__"                    
[113] "treedepth__"                    "n_leapfrog__"                  
[115] "divergent__"                    "energy__"                      

summary(fitbs)

 Family: student 
  Links: mu = identity; sigma = identity; nu = identity 
Formula: rate ~ (rtref * exp(e/bk * (1/tref - 1/(temp + 273.15))))/(1 + exp(eh/bk * (1/(th + 273.15) - 1/(temp + 273.15)))) 
         rtref ~ 1 + (1 | area)
         e ~ 1 + (1 | area)
         eh ~ 1 + (1 | area)
         th ~ 1 + (1 | area)
   Data: dat (Number of observations: 343) 
  Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 10000

Group-Level Effects: 
~area (Number of levels: 11) 
                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(rtref_Intercept)     0.14      0.16     0.00     0.58 1.00     2574     4177
sd(e_Intercept)         0.21      0.13     0.01     0.52 1.00     2850     3740
sd(eh_Intercept)        0.16      0.18     0.01     0.64 1.00     2425     2599
sd(th_Intercept)        0.48      0.57     0.02     2.02 1.00     2514     3850

Population-Level Effects: 
                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
rtref_Intercept     1.55      1.24     0.30     4.78 1.00     1878     2279
e_Intercept         1.89      0.93     0.36     3.71 1.00     1931     2368
eh_Intercept        2.21      0.83     0.67     3.89 1.00     2309     2601
th_Intercept        4.43      3.05     0.56    12.77 1.00     1863     2153

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     0.04      0.00     0.03     0.04 1.00    12112     8196
nu       10.37      4.30     5.10    21.61 1.00    11772     7771

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

pop_pars <- summary(fitbs)$fixed |> 
  as.data.frame() |> 
  rownames_to_column(var = "parameter") |> 
  mutate(parameter = str_remove(parameter, "_Intercept"),
         parameter = ifelse(parameter == "rtref", "r_tref", parameter)) |> 
  rename(pop_Estimate = Estimate) |> 
  dplyr::select(parameter, pop_Estimate)

mutate: changed 4 values (100%) of 'parameter' (0 new NA)

rename: renamed one variable (pop_Estimate)

pop_pars <- fitbs |> 
  gather_draws(b_rtref_Intercept,
               b_e_Intercept,
               b_eh_Intercept,
               b_th_Intercept) |> 
  rename(parameter = .variable) |> 
  mutate(parameter = str_remove(parameter, "b_"),
         parameter = str_remove(parameter, "_Intercept"),
         parameter = ifelse(parameter == "rtref", "r_tref", parameter)) |> 
  rename(.pop_value = .value )

rename: renamed one variable (parameter)

mutate (grouped): changed 40,000 values (100%) of 'parameter' (0 new NA)

rename: renamed one variable (.pop_value)

brms_pars <- fitbs |> 
  gather_draws(r_area__rtref[area, Intercept],
               r_area__e[area, Intercept],
               r_area__eh[area, Intercept],
               r_area__th[area, Intercept]) |> 
  rename(parameter = .variable) |> 
  mutate(parameter = str_remove(parameter, "r_area__"),
         parameter = ifelse(parameter == "rtref", "r_tref", parameter)) |> 
  dplyr::select(-Intercept) |>
  left_join(pop_pars, by = c(".chain", ".iteration", ".draw", "parameter")) |> 
  mutate(.value = .pop_value + .value) |> 
  summarise(Estimate = median(.value),
            lwr = quantile(.value, probs = 0.025),
            upr = quantile(.value, probs = 0.975)) |> 
  ungroup() |> 
  mutate(source = "brms")

rename: renamed one variable (parameter)

mutate (grouped): changed 440,000 values (100%) of 'parameter' (0 new NA)

Adding missing grouping variables: `Intercept`
left_join: added one column (.pop_value)
> rows only in x 0
> rows only in y ( 0)
> matched rows 440,000
> =========
> rows total 440,000
mutate (grouped): changed 440,000 values (100%) of '.value' (0 new NA)
summarise: now 44 rows and 6 columns, 2 group variables remaining (area,
Intercept)
ungroup: no grouping variables
mutate: new variable 'source' (character) with one unique value and 0% NA

pars <- bind_rows(brms_pars, nls_pars)

# Need to trim errorbars...
ggplot(pars, aes(area, Estimate, color = source)) + 
  geom_point() +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0) + 
  facet_wrap(~parameter, scales = "free")

pars |> 
  ggplot(aes(area, Estimate, color = source)) + 
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = lwr, ymax = upr), width = 0,
                position = position_dodge(width = 0.5)) + 
  coord_cartesian(ylim = c(-3, 20)) + 
  facet_wrap(~parameter, scales = "free")

# And without errorbars...
pars |> 
  ggplot(aes(area, Estimate, color = source)) + 
  geom_point(position = position_dodge(width = 0.5)) +
  facet_wrap(~parameter, scales = "free")

pars |> as.data.frame()

    area Intercept parameter     Estimate           lwr          upr source
1     BS Intercept         e 2.061649e+00  4.947591e-01 3.964354e+00   brms
2     BS Intercept        eh 2.139810e+00  3.325379e-01 3.878722e+00   brms
3     BS Intercept    r_tref 1.079088e+00  2.841505e-01 4.604358e+00   brms
4     BS Intercept        th 3.414204e+00  3.637865e-01 1.261392e+01   brms
5     BT Intercept         e 1.930650e+00  3.717316e-01 3.797633e+00   brms
6     BT Intercept        eh 2.108886e+00  6.432173e-01 3.838065e+00   brms
7     BT Intercept    r_tref 1.182941e+00  3.083056e-01 4.938161e+00   brms
8     BT Intercept        th 3.674797e+00  6.518049e-01 1.360176e+01   brms
9     FB Intercept         e 1.842479e+00  3.512824e-01 3.721149e+00   brms
10    FB Intercept        eh 2.151314e+00  6.014750e-01 3.861840e+00   brms
11    FB Intercept    r_tref 1.212869e+00  3.113947e-01 5.048426e+00   brms
12    FB Intercept        th 3.746232e+00  6.948678e-01 1.362966e+01   brms
13    FM Intercept         e 1.743456e+00  1.566586e-01 3.635704e+00   brms
14    FM Intercept        eh 2.189398e+00  7.287063e-01 3.918864e+00   brms
15    FM Intercept    r_tref 1.186454e+00  3.008972e-01 4.931285e+00   brms
16    FM Intercept        th 3.626702e+00  6.184890e-01 1.296625e+01   brms
17    HO Intercept         e 1.994180e+00  4.728469e-01 3.850492e+00   brms
18    HO Intercept        eh 2.211189e+00  3.330007e-01 3.957714e+00   brms
19    HO Intercept    r_tref 1.036137e+00  2.731101e-01 4.423051e+00   brms
20    HO Intercept        th 3.289899e+00  2.281909e-01 1.237871e+01   brms
21    JM Intercept         e 1.818157e+00  3.237866e-01 3.686721e+00   brms
22    JM Intercept        eh 2.182567e+00  6.983269e-01 3.874508e+00   brms
23    JM Intercept    r_tref 1.193769e+00  3.030681e-01 4.968800e+00   brms
24    JM Intercept        th 3.674315e+00  6.752544e-01 1.296660e+01   brms
25    MU Intercept         e 1.795793e+00  2.273802e-01 3.664792e+00   brms
26    MU Intercept        eh 2.212908e+00  5.996936e-01 3.939457e+00   brms
27    MU Intercept    r_tref 1.153410e+00  2.928628e-01 4.825046e+00   brms
28    MU Intercept        th 3.531675e+00  5.317666e-01 1.283821e+01   brms
29    RA Intercept         e 1.966905e+00  4.935831e-01 3.847171e+00   brms
30    RA Intercept        eh 2.117163e+00  4.394163e-01 3.829591e+00   brms
31    RA Intercept    r_tref 1.186353e+00  3.068444e-01 4.992200e+00   brms
32    RA Intercept        th 3.674079e+00  6.474994e-01 1.381782e+01   brms
33 SI_EK Intercept         e 1.778027e+00  1.495447e-01 3.669992e+00   brms
34 SI_EK Intercept        eh 2.262188e+00  5.736340e-01 3.992170e+00   brms
35 SI_EK Intercept    r_tref 1.082631e+00  2.758760e-01 4.560410e+00   brms
36 SI_EK Intercept        th 3.357972e+00  3.066126e-01 1.202503e+01   brms
37 SI_HA Intercept         e 1.717698e+00  3.254821e-02 3.588196e+00   brms
38 SI_HA Intercept        eh 2.214496e+00  7.216116e-01 3.951380e+00   brms
39 SI_HA Intercept    r_tref 1.156917e+00  2.916271e-01 4.812762e+00   brms
40 SI_HA Intercept        th 3.531649e+00  5.409756e-01 1.254292e+01   brms
41    TH Intercept         e 1.854305e+00  2.895028e-01 3.758030e+00   brms
42    TH Intercept        eh 2.159259e+00  5.607438e-01 3.883714e+00   brms
43    TH Intercept    r_tref 1.185886e+00  3.070300e-01 4.941389e+00   brms
44    TH Intercept        th 3.661620e+00  6.494001e-01 1.333225e+01   brms
45    JM      <NA>    r_tref 5.265193e-01 -1.770722e+00 2.823761e+00    nls
46    JM      <NA>         e 1.408110e+00 -3.406946e+00 6.223167e+00    nls
47    JM      <NA>        eh 2.375225e+00  8.382148e-01 3.912236e+00    nls
48    JM      <NA>        th 7.333556e+00 -1.463240e+01 2.929951e+01    nls
49    FM      <NA>    r_tref 2.513226e-01 -8.848557e-01 1.387501e+00    nls
50    FM      <NA>         e 9.504990e-02 -4.102377e+00 4.292477e+00    nls
51    FM      <NA>        eh 1.509774e+00 -1.200467e+01 1.502422e+01    nls
52    FM      <NA>        th 1.685407e+01 -6.856950e+01 1.022776e+02    nls
53 SI_HA      <NA>    r_tref 5.265193e-01 -2.179154e+02 2.189685e+02    nls
54 SI_HA      <NA>         e 0.000000e+00 -1.273175e+02 1.273175e+02    nls
55 SI_HA      <NA>        eh 5.164810e-01 -6.593563e+01 6.696859e+01    nls
56 SI_HA      <NA>        th 2.749949e+00 -8.122367e+03 8.127867e+03    nls
57 SI_EK      <NA>    r_tref 5.265193e-01 -8.040449e+00 9.093488e+00    nls
58 SI_EK      <NA>         e 1.795502e+00 -2.077437e+01 2.436538e+01    nls
59 SI_EK      <NA>        eh 3.079775e+00 -7.127788e+00 1.328734e+01    nls
60 SI_EK      <NA>        th 6.556893e+00 -4.312808e+01 5.624187e+01    nls
61    FB      <NA>    r_tref 5.265193e-01 -1.318479e+02 1.329010e+02    nls
62    FB      <NA>         e 6.408831e-01 -9.991975e+01 1.012015e+02    nls
63    FB      <NA>        eh 8.064473e-01 -1.525330e+01 1.686619e+01    nls
64    FB      <NA>        th 6.575958e+00 -3.844943e+03 3.858095e+03    nls
65    BT      <NA>    r_tref 2.428395e-01  1.290965e-01 3.565826e-01    nls
66    BT      <NA>         e 6.781416e-01 -9.498855e-01 2.306169e+00    nls
67    BT      <NA>        eh 3.083687e+00 -2.478754e+00 8.646128e+00    nls
68    BT      <NA>        th 1.551842e+01  5.400771e+00 2.563606e+01    nls
69    MU      <NA>    r_tref 2.540028e-01 -2.479580e+00 2.987585e+00    nls
70    MU      <NA>         e 4.668099e-01 -1.545223e+01 1.638585e+01    nls
71    MU      <NA>        eh 4.353853e+00 -1.042919e+02 1.129996e+02    nls
72    MU      <NA>        th 9.905223e+00 -1.712794e+01 3.693839e+01    nls
73    HO      <NA>    r_tref 3.487259e-01 -7.142516e+01 7.212261e+01    nls
74    HO      <NA>         e 4.871650e-01 -1.548278e+02 1.558022e+02    nls
75    HO      <NA>        eh 1.974064e-05 -3.106484e+02 3.106485e+02    nls
76    HO      <NA>        th 1.257663e+01 -1.843269e+08 1.843269e+08    nls
77    BS      <NA>    r_tref 1.787827e-01  1.569488e-01 2.006167e-01    nls
78    BS      <NA>         e 1.333788e+00  4.217272e-01 2.245849e+00    nls
79    BS      <NA>        eh 2.000000e+01 -3.126065e+01 7.126065e+01    nls
80    BS      <NA>        th 1.179472e+01  1.137974e+01 1.220969e+01    nls
81    RA      <NA>    r_tref 5.265193e-01 -1.953796e+05 1.953807e+05    nls
82    RA      <NA>         e 1.161895e+00 -1.847792e+04 1.848024e+04    nls
83    RA      <NA>        eh 1.012458e-01 -3.218889e+03 3.219092e+03    nls
84    RA      <NA>        th 1.685407e+01 -5.646609e+07 5.646613e+07    nls
85    TH      <NA>    r_tref 4.741047e-01 -3.312802e+03 3.313750e+03    nls
86    TH      <NA>         e 1.302458e-03 -2.433540e+03 2.433543e+03    nls
87    TH      <NA>        eh 1.494928e-04 -4.867156e+03 4.867157e+03    nls
88    TH      <NA>        th 1.685407e+01 -8.541512e+08 8.541512e+08    nls

Extra analysis

knitr::knit_exit()